# install.packages("remotes")
remotes::install_github("kadyb/rgugik")Skrypty geoprzetwarzania
Analizy geoprzestrzenne
library("terra")
library("rgugik")Główny Urząd Geodezji i Kartografii jest istotnym źródłem danych przestrzennych dla Polski. Dane można przeglądać i pobrać z Geoportalu lub wykorzystując różne usługi. W otwartych zbiorach danych znajdziemy m. in.:
- ortofotomapy,
- cyfrowe modele wysokościowe (CMW):
- numeryczny model terenu (NMT),
- numeryczny model pokrycia terenu (NMPT),
- chmury punktów,
- modele 3D budynków,
- Państwowy Rejestr Granic (PRG),
- Baza Danych Obiektów Topograficznych (BDOT),
- i inne.
Wyszukiwanie i pobieranie wymienionych zbiorów danych umożliwia pakiet rgugik.
Ortofotomapa
Ortofotomapa to rastrowe, ortogonalne i kartometryczne przedstawienie powierzchni terenu powstałe w wyniku cyfrowego przetwarzania zdjęć lotniczych lub satelitarnych. Podczas ortorektyfikacji usuwane są zniekształcenia geometryczne wynikające z rzeźby terenu przy użyciu cyfrowych modeli wysokości. Ortofotomapa posiada georeferencje, co pozwala na określenie współrzędnych geograficznych dla każdej komórki obrazu.
Cechy ortofotomapy:
- Rozdzielczość przestrzenna – związana z rozmiarem najmniejszego obiektu, który może zostać wykryty przez czujnik i jest określana przez rozmiar komórki obrazu (piksel). Im mniejsza komórka, tym więcej szczegółów reprezentuje. Zbyt duży rozmiar oznacza, że poszczególne obiekty na zdjęciu przestają być rozpoznawalne.
- Kompozycja – obrazy analogowe są przedstawione w odcieniach szarości, natomiast obrazy cyfrowe mogą składać się z naturalnych kolorów (RGB) lub bliskiej podczerwieni (NIR).
Wyszukiwanie
Do wyszukania arkuszy ortofotomapy służy funkcja ortho_request(), która jako argument przyjmuje geometrię. Jako przykładową geometrię możemy wykorzystać punkt. W tym celu należy stworzyć macierz, w której wierszy reprezentują punkty (w naszym przypadku tylko jeden), a kolumny współrzędne X i Y. Następnie należy dokonać konwersji do obiektu wektorowego używając funkcji vect() i definiując odpowiedni układ współrzędnych.
punkt = cbind(16.92, 52.41)
punkt = vect(punkt, crs = "EPSG:4326")dane = ortho_request(punkt)Możemy wyświetlić część otrzymanej ramki danych lub alternatywnie przeglądać całość używając funkcji View().
# wyświetl 10 pierwszych wierszy i 6 pierwszych kolumn
dane[1:10, 1:6] sheetID year resolution composition sensor CRS
1 N-33-130-D-d-1-2 2023 0.25 CIR Digital PL-1992
2 6.177.11.05.3.3 2023 0.05 CIR Digital PL-2000:S6
3 6.177.11.2 2023 0.50 RGB Digital PL-2000:S6
4 6.177.11.05.3.3 2023 0.05 RGB Digital PL-2000:S6
5 6.177.11.05.3 2023 0.10 CIR Digital PL-2000:S6
6 6.177.11.05.3 2023 0.10 RGB Digital PL-2000:S6
7 6.177.11.05 2023 0.20 CIR Digital PL-2000:S6
8 6.177.11.05 2023 0.20 RGB Digital PL-2000:S6
9 6.177.11.2 2023 0.50 CIR Digital PL-2000:S6
10 N-33-130-D-d-1-2 2021 0.25 CIR Digital PL-1992Standardowo dane możemy filtrować z uwzględnieniem zadanych parametrów.
dane[dane$year > 2016, ]
dane[dane$composition == "CIR", ]I sortować, np. według aktualności.
# kolejność malejąca (najnowsze dane)
dane[order(-dane$year), ]Pobieranie
Jako przykład pobierzmy dwie kompozycje tego samego obszaru wykonane w naturalnych barwach i z kanałem bliskiej podczerwieni z 2021 r. (ID: 75107_1047046_N-33-130-D-d-1-2 i 75106_1052150_N-33-130-D-d-1-2).
id = c("75107_1047046_N-33-130-D-d-1-2", "75106_1052150_N-33-130-D-d-1-2")
dane_sel = dane[dane$filename %in% id, ]Po selekcji potrzebnych danych, można je pobrać wykorzystując funkcję tile_download(). Możliwe jest również wskazanie katalogu, do którego powinny zostać pobrane obrazy (argument outdir).
Zazwyczaj warto zwiększyć domyślną wartość przekroczenia czasu połączenia (timeout) z domyślnych 60 sekund w przypadku dużych plików lub wolnego połączenia.
options(timeout = 600)
tile_download(dane_sel, outdir = "dane")Do wylistowania pobranych plików służy funkcja list.files(). Należy wskazać jakie pliki chcemy wczytać (pattern = "\\.tif$") i zapobiegawczo zwrócić pełne ścieżki do plików (full.names = TRUE).
pliki = list.files("dane", pattern = "\\.tif$", full.names = TRUE)
pliki[1] "dane/75106_1052150_N-33-130-D-d-1-2.tif"
[2] "dane/75107_1047046_N-33-130-D-d-1-2.tif"W ostatnim kroku możemy kolejno wczytać rastry i je wyświetlić.
# kompozycja w naturalnych barwach
r1 = rast(pliki[1])
plot(r1)
# kompozycja z bliską podczerwienią
r2 = rast(pliki[2])
plot(r2)
Cyfrowe modele wysokościowe
Cyfrowe modele wysokościowe to modele opisujące powierzchnię terenu. Powstają w wyniku obróbki zdjęć lotniczych, skanowania laserowego (LiDAR), pomiarów geodezyjnych czy interferometrii radarowej (InSAR). CMW są jednym z kluczowych zbiorów danych w systemach informacji geograficznej i stanowią podstawę wielu środowiskowych analiz przestrzennych. Ponadto są źródłem produktów pochodnych, takich jak cieniowanie, nachylenie czy szorstkość terenu.
CMW to ogólna nazwa grupy modeli o różnych cechach, uwzględniając:
- Numeryczny model terenu (Digital Terrain Model) – reprezentacja pozbawiona jakichkolwiek obiektów nad powierzchnią terenu, takich jak budynki czy drzewa.
- Numeryczny model pokrycia terenu (Digital Surface Model) – reprezentacja terenu wraz z jego pokryciem.
https://commons.wikimedia.org/w/index.php?title=File:DTM_DSM.svg
Cechy CMW:
- Format – można wyróżnić trzy główne struktury: GRID (regularna siatka punktów / komórek), TIN (nieregularna topologiczna sieć trójkątów) oraz linie konturowe (dane wektorowe). Obecnie najczęściej używanym formatem jest GRID.
- Dokładność – związana jest z pionowym błędem pomiaru.
- Rozdzielczość przestrzenna – związana jest z rozmiarem najmniejszego obiektu, który może zostać wykryty przez czujnik i jest określana przez rozmiar komórki obrazu (piksel). Im większa komórka, tym bardziej uogólniona forma terenu.
Meteoryt Morasko
Naszym obszarem analiz jest rezerwat przyrody Meteoryt Morasko położony w północnej części Poznania. Został on utworzony w 1976 roku w celu ochrony obszaru kraterów uderzeniowych, które według badaczy powstały w wyniku upadku meteorytu Morasko około 5 tysięcy lat temu. Ponadto ochroną objęty jest las grądowy z rzadkimi gatunkami roślin i ptaków.
https://naukawpolsce.pl/aktualnosci/news%2C402631%2Csto-lat-temu-odkryto-pierwszy-fragment-meteorytu-morasko.html
(Fot. PAP © 2012 / Jakub Kaczmarczyk)Więcej informacji znajdziesz tutaj:
Centroid (środek geometryczny) rezerwatu Morasko znajduje się na 16,895° długości geograficznej (X) i 52,489° szerokości geograficznej (Y).
wspolrzedne = matrix(c(16.895, 52.489), ncol = 2)
centroid = vect(wspolrzedne, type = "points", crs = "EPSG:4326")
centroid class : SpatVector
geometry : points
dimensions : 1, 0 (geometries, attributes)
extent : 16.895, 16.895, 52.489, 52.489 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) Dokonajmy również konwersji ze współrzędnych geograficznych na układ metryczny EPSG:2180.
centroid = project(centroid, "EPSG:2180")W kolejnym kroku stwórzmy przybliżoną strefę, która obejmie obszar rezerwatu.
bufor = buffer(centroid, width = 400)Stworzyliśmy bufor o szerokości 400 m, a teraz przygotujmy prostą wizualizację.
plot(bufor, main = "Bufor rezerwatu Morasko")
plot(centroid, col = "blue", add = TRUE)
Następnie możemy wyszukać dostępne dane wysokościowe dla tego obszaru za pomocą funkcji DEM_request() (jest ona analogiczna do funkcji ortho_request()). Jako argument musimy wskazać nasz utworzony poligon.
dane = DEM_request(centroid)Oczywiście pozyskane dane możemy sprawdzić wywołując obiekt dane lub przejrzeć je interaktywnie używając funkcji View().
# wyświetl 10 pierwszych wierszy i 6 pierwszych kolumn
dane[1:10, 1:6] sheetID year product format source
1 6.179.11.14 2021 DSM ARC/INFO ASCII GRID Laser scanning
2 6.179.11.14 2021 DTM ARC/INFO ASCII GRID Laser scanning
3 N-33-130-D-b-1-1 2021 DTM ARC/INFO ASCII GRID Aerial photo
4 N-33-130-D-b-1-1 2020 DTM ARC/INFO ASCII GRID Aerial photo
5 N-33-130-D-b-1-1-4-1 2019 PointCloud LAS Laser scanning
6 N-33-130-D-b-1-1 2012 DSM ASCII XYZ GRID Laser scanning
7 N-33-130-D-b-1-1 2012 DSM ARC/INFO ASCII GRID Laser scanning
8 N-33-130-D-b-1-1 2012 DTM ASCII XYZ GRID Laser scanning
9 N-33-130-D-b-1-1-4-1 2012 PointCloud LAS Laser scanning
10 N-33-130-D-b-1-1 2012 DTM ARC/INFO ASCII GRID Laser scanning
resolution
1 0.50 m
2 0.50 m
3 5.00 m
4 5.00 m
5 12 p/m2
6 0.50 m
7 0.50 m
8 1.00 m
9 12 p/m2
10 1.00 mJak możemy zauważyć powyższe metadane opisują produkty o różnych formatach, aktualności, rozdzielczości oraz dokładności. Do naszej analizy potrzebujemy numerycznego modelu terenu (DTM) i numerycznego modelu pokrycia terenu (DSM) w formacie “ARC/INFO ASCII GRID”. Dokonajmy selekcji danych, tworząc dwie ramki danych i następnie łącząc je ze sobą przy pomocy funkcji rbind().
DTM_sel = dane[dane$format == "ARC/INFO ASCII GRID" &
dane$product == "DTM" &
dane$year == 2019, ]
DSM_sel = dane[dane$format == "ARC/INFO ASCII GRID" &
dane$product == "DSM" &
dane$year == 2019, ]
# połączenie powyższych ramek danych
dane_sel = rbind(DTM_sel, DSM_sel)
dane_sel[, 1:6] sheetID year product format source resolution
19 N-33-130-D-b-1-1 2019 DTM ARC/INFO ASCII GRID Laser scanning 1.00 m
15 N-33-130-D-b-1-1 2019 DSM ARC/INFO ASCII GRID Laser scanning 0.50 mWykorzystajmy funkcję tile_download() do pobrania tych dwóch produktów.
options(timeout = 600)
tile_download(dane_sel, outdir = "dane")Przetwarzanie
Po pobraniu danych możemy je wczytać do sesji.
DTM = rast("dane/73044_917579_N-33-130-D-b-1-1.asc")
DSM = rast("dane/73043_917495_N-33-130-D-b-1-1.asc")Możemy również nadać warstwom odpowiednie nazwy oraz, co ważniejsze, przypisać im układy współrzędnych zdefiniowane w obiekcie dane_sel (atrybut CRS).
# nadanie nazw warstwom
names(DTM) = "DTM"
names(DSM) = "DSM"
# ustawienie układu współrzędnych
crs(DTM) = crs(DSM) = "EPSG:2180"Podczas pobierania prawdopodobnie zauważyłeś, że rastry różnią się czterokrotnie rozmiarem. Wynika to z różnicy ich rozdzielczości przestrzennej. W takiej sytuacji należy ujednolicić je do jednakowej rozdzielczości, aby móc je połączyć (nałożyć na siebie). Znacznie lepiej jest użyć niższej rozdzielczości niż ją sztucznie zwiększać, ponieważ nie możemy uzyskać więcej informacji, a przetwarzanie będzie szybsze. W tym celu użyjmy funkcji resample() i zapiszmy wynik na dysku określając ścieżkę w argumencie filename.
DSM = resample(DSM, DTM, method = "near", filename = "dane/DSM_1.tif")Teraz oba modele mają te same wymiary (liczbę wierszy i kolumn) oraz rozdzielczość przestrzenną. Możemy więc połączyć je w jeden obiekt o nazwie DEM.
DEM = c(DTM, DSM)
nlyr(DEM)[1] 2Używając funkcji crop() możemy ograniczyć zasięg przestrzenny (bounding box) do zasięgu przestrzennego rezerwatu Morasko. W celu wykluczenia z analizy wartości poza obszarem poligonu, należy zastosować funkcję mask(). Można tę czynność również wykonać w jednym kroku, używają funkcji crop() z argumentem mask = TRUE.
DEM_crop = crop(DEM, bufor)
DEM_mask = mask(DEM_crop, bufor)Zauważ, że domyślnie kolory zielony reprezentuje wysokie wartości, a kolor pomarańczowy niskie. Do wizualizacji danych topograficznych warto odwrócić paletę kolorów, tj. terrain.colors(99, alpha = NULL).
par(mfrow = c(1, 2)) # wyświetl 2 rastry obok siebie
plot(DEM_crop$DTM, main = "Docięcie", col = terrain.colors(99, alpha = NULL))
plot(bufor, add = TRUE)
plot(DEM_mask$DTM, main = "Maskowanie", col = terrain.colors(99, alpha = NULL))
plot(bufor, add = TRUE)
W pierwszej ćwiartce okręgu widzimy pięć mniejszych okręgów. Są to kratery powstałe po uderzeniu meteorytu Morasko. Największy znaleziony fragment waży 272 kg i jest największym meteorytem znalezionym w Polsce. Kolekcję można zobaczyć w Muzeum Ziemi WNGIG w Poznaniu.
Obliczmy szerokość krateru na podstawie poprzecznego profilu terenu. Przyjmijmy, że punkt A ma współrzędne (357280, 515980), a punkt B (357122, 515760).
punkty = matrix(c(357280, 515980,
357122, 515760),
ncol = 2, byrow = TRUE)
linia = vect(punkty, type = "lines", crs = "EPSG:2180")
linia class : SpatVector
geometry : lines
dimensions : 1, 0 (geometries, attributes)
extent : 357122, 357280, 515760, 515980 (xmin, xmax, ymin, ymax)
coord. ref. : ETRF2000-PL / CS92 (EPSG:2180) plot(DEM_mask$DTM, main = "DTM [m]", col = terrain.colors(99, alpha = NULL))
plot(linia, col = "red", lwd = 2, add = TRUE)
text("A", x = 357320, y = 515980, cex = 0.8)
text("B", x = 357100, y = 515760, cex = 0.8)
W kolejnym kroku pozyskamy wartości wysokości dla wyznaczonego profilu za pomocą funkcji extract().
profil = extract(DEM, linia)
profil = profil$DTM
plot(profil, type = "l", xlab = "Indeks komórki", ylab = "Wysokość n.p.m. [m]")
Na profilu widoczny jest szum wynikający z metody pozyskania danych. Można to wygładzić w prosty sposób wykorzystując funkcję loess().
profil = loess(profil ~ seq_along(profil), span = 0.1)
profil = profil$fitted
plot(profil, type = "l", xlab = "Indeks komórki", ylab = "Wysokość n.p.m. [m]")
Powyższe ryciny na osi X zawierają indeks komórki. Lepszym rozwiązaniem będzie przedstawienie tej osi jako odległości od punktu startowego (A). Aby to wykonać musimy obliczyć średnią odległość między komórkami, czyli podzielić długość linii przez liczbę komórek.
odleglosc = perim(linia) / length(profil)
odleglosc[1] 0.7146646Następnie stwórzmy wektor odległości kolejnych punktów od punktu startowego wykorzystując sumę skumulowaną cumsum().
odleglosc = cumsum(rep(odleglosc, length(profil)))
odleglosc[1:5] # odległość pierwszych 5 punktów[1] 0.7146646 1.4293293 2.1439939 2.8586585 3.5733232Po tej operacji, możemy zaprezentować finalną wersję wykresu z zaznaczonym zasięgiem największego krateru uderzeniowego wynoszącą około 90 m.
plot(profil ~ odleglosc, type = "l", xlab = "Odległość [m]",
ylab = "Wysokość n.p.m. [m]", main = "Numeryczny model terenu")
abline(v = c(50, 140), col = "blue")
W ostatnim kroku sprawdźmy wysokość obiektów (na tym obszarze są to drzewa). W tym celu należy odjąć NMT od NMPT. Różnica nazywana jest znormalizowanym NMPT, ponieważ przyjmuje wysokość terenu jako odniesienie.
nDSM = DEM$DSM - DEM$DTM
nDSMclass : SpatRaster
dimensions : 2379, 2188, 1 (nrow, ncol, nlyr)
resolution : 1, 1 (x, y)
extent : 355732.5, 357920.5, 514649.5, 517028.5 (xmin, xmax, ymin, ymax)
coord. ref. : ETRF2000-PL / CS92 (EPSG:2180)
source(s) : memory
varname : 73044_917579_N-33-130-D-b-1-1
name : DSM
min value : -4.159996
max value : 67.769989 # nadpisz wartości poniżej 0
nDSM[nDSM < 0] = 0Powyższą operację odejmowania można wykonać również używając funkcji lapp().
plot(nDSM, main = "Wysokość drzew [m]",
col = hcl.colors(9, palette = "Greens", rev = TRUE))
Zadanie
10. Pobierz minimum dwa sąsiadujące ze sobą kafelki ortofotomapy z tej samej serii i połącz je:
- do jednego pliku .tiff używając funkcji
merge(), - do jednego wirtualnego pliku .vrt używając funkcji
vrt().
Sprawdź zajmowaną ilość miejsca przez te pliki na dysku wykorzystując funkcję file.size() (wynik zwracany jest w bajtach). Sprawdź również zawartość pliku .vrt (czym on jest w rzeczywistości?). Następnie, zmniejsz rozdzielczość mozaiki do 10 m i zapisz wynik. Jak zmieniła się jakość w stosunku do oryginalnego zdjęcia?
11. Wykonaj analizę ukształtowania terenu wybranego obszaru, w której uwzględnisz:
- wizualizacje NMT i NMPT,
- statystyki opisowe wysokości terenu (wartość minimalna, maksymalna, średnia oraz odchylenie standardowe),
- profil wysokościowy,
- wizualizacje wysokości obiektów (znormalizowany NMPT).
Pamiętaj, że:
- Rastry powinny posiadać taką samą rozdzielczość przestrzenną oraz układ współrzędnych.
- Różnica między NMPT a NMT nie może być ujemna.